Mini-Project #03: Visualizing and Maintaining the Green Canopy of NYC

Author

Hajara Muzammal

Introduction

New York City’s Department of Parks and Recreation manages nearly 900,000 trees across all five boroughs, providing critical environmental benefits to millions of residents, yet access to this green infrastructure remains unequal across neighborhoods. This project analyzes the NYC TreeMap dataset alongside City Council District boundaries to explore patterns in tree distribution, health, and species diversity at the district level. Using spatial analysis and data visualization, we identify opportunities for targeted intervention and propose a data-driven tree revitalization program that addresses gaps in the urban canopy and promotes environmental equity.

Data Preparation

Show code
# Load required packages
library(tidyverse)
library(sf)
library(httr2)
#| message: false
#| warning: false
#| output: false

Task 1: Download NYC City Council District Boundaries

Show code
suppressPackageStartupMessages({
  library(sf)
  library(fs)
})

NYC_City_Council <- function(url){
  mp03 <- file.path("data", "mp03")
  if (!dir.exists(mp03)) {
    dir.create(mp03, showWarnings=FALSE, recursive=TRUE)
  }
  
  zip_path <- file.path(mp03, "NYC City Council District Boundaries.zip")
  if (!file.exists(zip_path)) {
    download.file(url,
                  destfile = zip_path,
                  mode = "wb")
  }
  
  shp_file <- dir_ls(mp03, recurse = TRUE, glob = "*.shp")
  if (length(shp_file) == 0) {
    unzip(zip_path, exdir = mp03)
    shp_file <- dir_ls(mp03, recurse = TRUE, glob = "*.shp")
  }
  
  NYC_file <- st_read(shp_file[1], quiet = TRUE)
  NYC_file <- st_transform(NYC_file, crs = "WGS84")
  
  return(NYC_file)
}

council <- NYC_City_Council("https://s-media.nyc.gov/agencies/dcp/assets/files/zip/data-tools/bytes/city-council/nycc_25c.zip")

Task 2: Download Tree Points

Show code
#| label: download-tree-data
#| include: false
#| cache: true

library(httr2)
library(sf)
library(dplyr)

# Function to download NYC Tree Points using httr2
download_tree_points <- function(
  base_url = "https://data.cityofnewyork.us/resource/hn5i-inap.geojson",
  limit = 5000,
  data_dir = "data/mp03/trees",
  pause = 1
) {
  
  # Create trees directory if it doesn't exist
  if (!dir.exists(data_dir)) {
    dir.create(data_dir, recursive = TRUE)
    message("Created directory: ", data_dir)
  }
  
  files <- character()
  offset <- 0L
  page <- 1L
  
  message("Downloading NYC Tree Points data...")
  
  repeat {
    # Construct filename with consistent naming schema
    fname <- file.path(data_dir, sprintf("treepoints_%08d_%05d.geojson", offset, limit))
    
    # Only download if file doesn't exist
    if (!file.exists(fname)) {
      message("Downloading page ", page, " (offset: ", offset, ")...")
      
      # Build request using httr2
      req <- request(base_url) |>
        req_url_query(
          `$limit`  = limit,
          `$offset` = offset
        ) |>
        req_timeout(300) |>
        req_retry(max_tries = 3)
      
      # Perform request
      resp <- req_perform(req)
      writeBin(resp_body_raw(resp), fname)
      message("  Saved to: ", basename(fname))
      Sys.sleep(pause)  # Be respectful to the API
      
    } else {
      message("Page ", page, " already exists, skipping")
    }
    
    # Read the file to check how many records we got
    trees_page <- try(st_read(fname, quiet = TRUE), silent = TRUE)
    nret <- if (inherits(trees_page, "sf")) nrow(trees_page) else 0L
    
    message("  Page ", page, " contains ", nret, " records")
    files <- c(files, fname)
    
    # If we got fewer records than limit, we've reached the end
    if (nret < limit) {
      message("\nReached end of dataset at page ", page)
      break
    }
    
    offset <- offset + limit
    page <- page + 1L
  }
  
  # Read all GeoJSON files using st_read
  message("\nReading and combining ", length(files), " GeoJSON files...")
  
  all_trees <- lapply(files, function(f) {
    df <- st_read(f, quiet = TRUE)
    
    # Convert all columns to character except geometry
    geom_col <- attr(df, "sf_column")
    char_cols <- setdiff(names(df), geom_col)
    
    for (col in char_cols) {
      df[[col]] <- as.character(df[[col]])
    }
    
    return(df)
  })
  
  # Combine all datasets using bind_rows
  trees_combined <- bind_rows(all_trees)
  
  # Convert columns to proper types
  if ("planteddate" %in% names(trees_combined)) {
    trees_combined$planteddate <- as.POSIXct(
      trees_combined$planteddate,
      format = "%Y-%m-%dT%H:%M:%S",
      tz = "UTC"
    )
  }
  
  numeric_cols <- c("tree_dbh", "dbh", "stump_diam", "zipcode", 
                    "latitude", "longitude", "council_district", "census_tract")
  for (col in numeric_cols) {
    if (col %in% names(trees_combined)) {
      trees_combined[[col]] <- as.numeric(trees_combined[[col]])
    }
  }
  
  message("Download and combination complete!")
  message("Total records: ", format(nrow(trees_combined), big.mark = ","))
  
  return(trees_combined)
}

# Download the tree data (or read from cache)
trees <- download_tree_points(limit = 5000, pause = 1)

Task 3: Plot All Tree Points

Data Integration and Initial Exploration

Show code
#| label: load-data
#| message: false
#| warning: false

library(tidyverse)
library(sf)
library(httr2)

# Load council districts from already downloaded files
council_districts <- st_read("data/mp03/geo_export_0231ea95-571c-48c6-869c-a4f70a7fa223.shp", quiet = TRUE) |>
  st_transform(crs = "WGS84")

# Load trees from already downloaded files
tree_files <- list.files("data/mp03/trees", pattern = "\\.geojson$", full.names = TRUE)

all_trees <- lapply(tree_files, function(f) {
  df <- st_read(f, quiet = TRUE)
  geom_col <- attr(df, "sf_column")
  char_cols <- setdiff(names(df), geom_col)
  for (col in char_cols) df[[col]] <- as.character(df[[col]])
  return(df)
})

trees <- bind_rows(all_trees)

# Convert to proper data types
if ("planteddate" %in% names(trees)) {
  trees$planteddate <- as.POSIXct(trees$planteddate, 
                                   format = "%Y-%m-%dT%H:%M:%S", tz = "UTC")
}

numeric_cols <- c("dbh", "zipcode", "latitude", "longitude")
for (col in numeric_cols) {
  if (col %in% names(trees)) {
    trees[[col]] <- as.numeric(trees[[col]])
  }
}

# Create spatial join
trees_with_districts <- st_join(trees, council_districts, join = st_intersects)

cat("Loaded:", nrow(council_districts), "council districts\n")
Loaded: 51 council districts
Show code
cat("Loaded:", format(nrow(trees), big.mark = ","), "trees\n")
Loaded: 1,934,587 trees
Show code
## TASK 3: PLOT ALL TREE POINTS 
#| label: task3-plot-trees
#| message: false
#| warning: false
#| fig-width: 10
#| fig-height: 8

library(ggplot2)
library(dplyr)
library(sf)

# Simplify council districts for faster plotting
council_simplified <- council_districts %>% 
  st_transform(crs = "WGS84") %>%
  mutate(geometry = st_simplify(geometry, dTolerance = 5))

# Transform trees to WGS84
trees_wgs84 <- st_transform(trees, crs = "WGS84")

# Sample trees for performance (use 50,000 for better visualization)
set.seed(123)
trees_sample <- trees_wgs84 %>% 
  slice_sample(n = 50000)

# Create the plot
ggplot() +
  geom_sf(data = council_simplified, 
          fill = "white", 
          color = "gray20", 
          linewidth = 0.3) +
  geom_sf(data = trees_sample,
          color = "forestgreen", 
          alpha = 0.3, 
          size = 0.5) +
  labs(
    title = paste("NYC Tree Map (", format(nrow(trees_sample), big.mark=","), " Trees)"),
    subtitle = "Overlay on City Council Districts",
    caption = "Data: NYC OpenData & NYC Planning"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    plot.subtitle = element_text(size = 12),
    panel.background = element_rect(fill = "white"),
    panel.grid = element_line(color = "gray95"),
    axis.text = element_blank(),
    axis.title = element_blank()
  )

Above is a visual of the NYC Council Districts.

Extra credit 1 - Visualization

Show code
# First, check if have the required objects
cat("Checking required objects:\n")
Checking required objects:
Show code
cat("trees exists:", exists("trees"), "\n")
trees exists: TRUE 
Show code
cat("council_districts exists:", exists("council_districts"), "\n")
council_districts exists: TRUE 
Show code
# If trees doesn't exist, reload it
if (!exists("trees")) {
  trees <- download_tree_points(limit = 5000, pause = 1)
}

# If council_districts doesn't exist, reload it
if (!exists("council_districts")) {
  council_districts <- download_council_districts()
}

# Now run the leaflet code
library(leaflet)
library(htmlwidgets)
library(sf)

# Create interactive map using leaflet
# Sample trees for performance (10,000 points)
set.seed(123)

trees_sample <- trees %>%
  slice_sample(n = 10000) %>%
  mutate(
    species_info = paste0("<b>", genusspecies, "</b><br>",
                          "Condition: ", tpcondition)
  )

# Extract coordinates for leaflet
tree_coords <- st_coordinates(trees_sample)
trees_sample$lon <- tree_coords[, "X"]
trees_sample$lat <- tree_coords[, "Y"]

# Create leaflet map
interactive_map <- leaflet(width = "100%", height = "600px") %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  
  # Add council district boundaries
  addPolygons(
    data = council_districts %>% st_transform(crs = 4326),
    fillColor = "transparent",
    color = "black",
    weight = 2,
    opacity = 0.7,
    label = ~paste("Council District", coun_dist),
    group = "Districts"
  ) %>%
  
  # Add tree points
  addCircleMarkers(
    data = trees_sample,
    lng = ~lon,
    lat = ~lat,
    radius = 3,
    color = "darkgreen",
    fillColor = "forestgreen",
    fillOpacity = 0.6,
    stroke = FALSE,
    popup = ~species_info,
    label = ~genusspecies,
    group = "Trees",
    clusterOptions = markerClusterOptions()
  ) %>%
  
  # Add layer controls
  addLayersControl(
    overlayGroups = c("Districts", "Trees"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  
  # Add scale bar
  addScaleBar(position = "bottomleft") %>%
  
  # Set initial view to NYC
  setView(lng = -73.935242, lat = 40.730610, zoom = 11)

interactive_map

Use your mouse wheel or pinch gesture to zoom, and click-drag to pan across the map. Click any tree marker to view detailed information including species, condition, and location, or simply hover over a marker for a quick species preview. The layer control panel in the upper-right corner allows you to toggle district boundaries and tree markers on or off. For improved performance, trees are automatically clustered at higher zoom levels—zoom in to view individual trees. Double-click anywhere on the map to reset to the default NYC view. The scale bar in the bottom-left helps estimate distances.

Task 4: District-Level Analysis of Tree Coverage

  1. Which council district has the most trees?
Show code
library(knitr)

district_tree_counts <- trees_with_districts %>%
  st_drop_geometry() %>%
  filter(!is.na(coun_dist)) %>%  # Changed to lowercase
  count(coun_dist, name = "tree_count") %>%
  arrange(desc(tree_count))

# Create formatted table for top 10 with row numbers
district_tree_counts %>%
  head(10) %>%
  mutate(Rank = row_number()) %>%
  select(Rank, coun_dist, tree_count) %>%  # Changed to lowercase
  kable(
    col.names = c("Rank", "Council District", "Number of Trees"),
    format.args = list(big.mark = ","),
    caption = "Table 1: Top 10 NYC Council Districts by Tree Count",
    align = c("c", "c", "c")
  )
Table 1: Top 10 NYC Council Districts by Tree Count
Rank Council District Number of Trees
1 51 130,877
2 19 90,448
3 50 90,118
4 23 81,337
5 13 60,441
6 49 60,066
7 31 57,174
8 32 53,470
9 27 51,736
10 39 51,116

District 51 has the most trees.

2. Which district has the highest density of trees?

Show code
library(knitr)

district_density <- trees_with_districts %>%
  st_drop_geometry() %>%
  filter(!is.na(coun_dist)) %>%  # Changed to lowercase
  count(coun_dist, name = "tree_count") %>%
  left_join(council_districts %>%  # Changed from 'council'
              st_drop_geometry() %>% 
              select(coun_dist, shape_area),  # Changed to lowercase
            by = "coun_dist") %>%  # Changed to lowercase
  mutate(tree_density = tree_count / shape_area) %>%  # Changed to lowercase
  arrange(desc(tree_density))

# Formatted table
district_density %>%
  head(10) %>%
  mutate(Rank = row_number()) %>%
  select(Rank, coun_dist, tree_count, tree_density) %>%  # Changed to lowercase
  kable(
    col.names = c("Rank", "Council District", "Number of Trees", "Tree Density"),
    format.args = list(big.mark = ","),
    caption = "Table 2: Top 10 NYC Council Districts by Tree Density",
    align = c("c", "c", "c", "c"),
    digits = c(0, 0, 0, 6)
  )
Table 2: Top 10 NYC Council Districts by Tree Density
Rank Council District Number of Trees Tree Density
1 39 51,116 0.000432
2 9 24,204 0.000392
3 2 20,926 0.000367
4 16 22,771 0.000366
5 14 19,104 0.000362
6 35 27,485 0.000346
7 7 28,632 0.000329
8 36 24,829 0.000326
9 44 31,893 0.000322
10 41 25,387 0.000320

District 39 has the highest tree density.

  1. Which district has highest fraction of dead trees?
Show code
library(knitr)

dead_tree_fraction <- trees_with_districts %>%
  st_drop_geometry() %>%
  filter(!is.na(coun_dist)) %>%  # Changed to lowercase
  group_by(coun_dist) %>%  # Changed to lowercase
  summarise(
    total_trees = n(),
    # Use tpstructure instead of status - "Retired" and "Stump" indicate dead trees
    dead_trees = sum(tpstructure %in% c("Retired", "Stump"), na.rm = TRUE),
    dead_fraction = dead_trees / total_trees,
    .groups = "drop"
  ) %>%
  arrange(desc(dead_fraction))

# Formatted table
dead_tree_fraction %>%
  head(10) %>%
  mutate(
    Rank = row_number(),
    dead_percentage = dead_fraction * 100
  ) %>%
  select(Rank, coun_dist, total_trees, dead_trees, dead_percentage) %>%  # Changed to lowercase
  kable(
    col.names = c("Rank", "Council District", "Total Trees", "Dead Trees", "Dead Tree %"),
    format.args = list(big.mark = ","),
    caption = "Table 3: Top 10 NYC Council Districts by Dead Tree Percentage",
    align = c("c", "c", "c", "c", "c"),
    digits = c(0, 0, 0, 0, 2)
  )
Table 3: Top 10 NYC Council Districts by Dead Tree Percentage
Rank Council District Total Trees Dead Trees Dead Tree %
1 32 53,470 13,211 24.71
2 2 20,926 5,033 24.05
3 48 31,194 7,492 24.02
4 1 22,004 5,217 23.71
5 3 21,310 4,720 22.15
6 11 43,713 9,522 21.78
7 4 20,441 4,441 21.73
8 30 43,282 9,367 21.64
9 14 19,104 4,115 21.54
10 16 22,771 4,893 21.49

District 16 has the highest fraction of dead trees.

  1. What is the most common tree species in Manhattan?
Show code
library(knitr)

# Add borough information based on district ranges
trees_with_borough <- trees_with_districts %>%
  mutate(Borough = case_when(
    coun_dist >= 1 & coun_dist <= 10 ~ "Manhattan",
    coun_dist >= 11 & coun_dist <= 18 ~ "Bronx",
    coun_dist >= 19 & coun_dist <= 32 ~ "Queens",
    coun_dist >= 33 & coun_dist <= 48 ~ "Brooklyn",
    coun_dist >= 49 & coun_dist <= 51 ~ "Staten Island",
    TRUE ~ NA_character_
  ))

manhattan_species <- trees_with_borough %>%
  st_drop_geometry() %>%
  filter(Borough == "Manhattan", !is.na(genusspecies)) %>%
  count(genusspecies, sort = TRUE)

# Formatted table
manhattan_species %>%
  head(10) %>%
  mutate(Rank = row_number()) %>%
  select(Rank, genusspecies, n) %>%
  kable(
    col.names = c("Rank", "Tree Species", "Count"),
    format.args = list(big.mark = ","),
    caption = "Table 4: Top 10 Most Common Tree Species in Manhattan",
    align = c("c", "l", "c")
  )
Table 4: Top 10 Most Common Tree Species in Manhattan
Rank Tree Species Count
1 Gleditsia triacanthos var. inermis - Thornless honeylocust 33,333
2 Platanus x acerifolia - London planetree 22,419
3 Pyrus calleryana - Callery pear 17,238
4 Quercus palustris - pin oak 15,292
5 Ginkgo biloba - maidenhair tree 14,339
6 Zelkova serrata - Japanese zelkova 10,690
7 Styphnolobium japonicum - Japanese pagoda tree 10,515
8 Tilia cordata - littleleaf linden 8,535
9 Unknown - Unknown 6,903
10 Ulmus americana - American elm 6,421

The most common tree species is honeylocust.

  1. What is the species of the tree closest to Baruch’s campus?
Show code
#| label: baruch-closest-tree
#| message: false
#| warning: false

library(knitr)
library(sf)
library(dplyr)

# Create helper function for point creation
new_st_point <- function(lon, lat) {
  st_sfc(st_point(c(lon, lat))) %>%
    st_set_crs("WGS84")
}

# Baruch College coordinates (Lexington Ave & E 24th St)
baruch_point <- new_st_point(-73.9838, 40.7402)

# Make sure trees are in WGS84
trees_wgs84 <- st_transform(trees, crs = "WGS84")

# Calculate distance from each tree to Baruch
trees_with_distance <- trees_wgs84 %>%
  mutate(distance = st_distance(geometry, baruch_point))

# Show TOP 5 closest trees
top_5_closest <- trees_with_distance %>%
  arrange(distance) %>%
  slice(1:5) %>%
  st_drop_geometry() %>%
  select(genusspecies, distance) %>%  # Changed from spc_common
  mutate(distance_m = round(as.numeric(distance), 1)) %>%
  select(genusspecies, distance_m)

top_5_closest %>%
  kable(
    col.names = c("Tree Species", "Distance (meters)"),
    caption = "Top 5 Trees Closest to Baruch College Campus",
    align = c("l", "c")
  )
Top 5 Trees Closest to Baruch College Campus
Tree Species Distance (meters)
Quercus acutissima - sawtooth oak 16.4
Quercus imbricaria - shingle oak 28.7
Quercus acutissima - sawtooth oak 28.7
Quercus acutissima - sawtooth oak 30.4
Quercus acutissima - sawtooth oak 30.4

Sawtooth oak is the closest tree closest to Baruch.

Task 5: NYC Parks Proposal

I will be analyzing all the data that is needed for this proposal.

Show code
# ============================================================================
# Task 5: NYC Parks Proposal - District Tree Program
# ============================================================================

library(tidyverse)
library(sf)
library(scales)

# ============================================================================
# STEP 1: Choose Your District
# ============================================================================

# Choose one of these options:
# Option 1: Baruch's district (District 2 - Manhattan)
my_district <- 2

# Option 2: Your district (replace with your district number)
# my_district <- XX

# Get your district data
my_district_boundary <- council_districts |>
  filter(coun_dist == my_district)

# Get trees in your district
my_district_trees <- trees_with_districts |>
  filter(coun_dist == my_district)
Show code
# ============================================================================
# STEP 2: Choose Your Project Focus
# ============================================================================

# Example: Replace Dead/Retired Trees and Stumps
# Count trees by structure type
tree_structure_summary <- my_district_trees |>
  st_drop_geometry() |>
  count(tpstructure, sort = TRUE)

print(tree_structure_summary)
       tpstructure     n
1             Full 15861
2          Retired  4453
3            Stump   580
4            Shaft    27
5 Stump - Uprooted     5
Show code
# Count stumps and retired trees
stumps_count <- my_district_trees |>
  st_drop_geometry() |>
  filter(tpstructure %in% c("Stump", "Retired")) |>
  nrow()

cat("\nStumps and retired trees in your district:", stumps_count, "\n")

Stumps and retired trees in your district: 5033 
Show code
# ============================================================================
# STEP 3: Quantitative Scope Statement
# ============================================================================

# Calculate project scope
project_scope <- my_district_trees |>
  st_drop_geometry() |>
  summarise(
    total_trees = n(),
    stumps = sum(tpstructure == "Stump", na.rm = TRUE),
    retired = sum(tpstructure == "Retired", na.rm = TRUE),
    poor_condition = sum(tpcondition == "Poor", na.rm = TRUE),
    to_replace = stumps + retired
  )

cat("\nProject Scope:\n")

Project Scope:
Show code
print(project_scope)
  total_trees stumps retired poor_condition to_replace
1       20926    580    4453            546       5033
Show code
# ============================================================================
# STEP 4: Zoomed-in Map of Your District
# ============================================================================

# Map showing all trees in your district
map_my_district <- ggplot() +
  # District boundary
  geom_sf(data = my_district_boundary, 
          fill = "lightgray", 
          color = "black", 
          linewidth = 1.5) +
  
  # All trees
  geom_sf(data = my_district_trees,
          aes(color = tpstructure),
          size = 0.5,
          alpha = 0.6) +
  
  scale_color_manual(
    name = "Tree Status",
    values = c("Full" = "darkgreen", 
               "Retired" = "orange", 
               "Stump" = "red")
  ) +
  
  theme_minimal() +
  labs(
    title = paste("Trees in Council District", my_district),
    subtitle = paste(format(nrow(my_district_trees), big.mark = ","), "total trees")
  ) +
  theme(
    axis.text = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    legend.position = "bottom"
  )

print(map_my_district)

Show code
# Alternative: Map showing only stumps and retired trees
map_replacement_targets <- ggplot() +
  geom_sf(data = my_district_boundary, 
          fill = "white", 
          color = "black", 
          linewidth = 1.5) +
  
  # Only stumps and retired trees
  geom_sf(data = my_district_trees |> 
            filter(tpstructure %in% c("Stump", "Retired")),
          aes(color = tpstructure),
          size = 1,
          alpha = 0.7) +
  
  scale_color_manual(
    name = "Replacement Needed",
    values = c("Retired" = "orange", "Stump" = "red")
  ) +
  
  theme_minimal() +
  labs(
    title = paste("Trees to Replace in District", my_district),
    subtitle = paste(stumps_count, "stumps and retired trees")
  ) +
  theme(
    axis.text = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    legend.position = "bottom"
  )

print(map_replacement_targets)

Show code
# ============================================================================
# STEP 5: Compare to Other Districts
# ============================================================================

# Compare stumps/retired trees across districts
district_comparison <- trees_with_districts |>
  st_drop_geometry() |>
  filter(!is.na(coun_dist)) |>
  group_by(coun_dist) |>
  summarise(
    total_trees = n(),
    stumps_retired = sum(tpstructure %in% c("Stump", "Retired"), na.rm = TRUE),
    percent_need_replacement = (stumps_retired / total_trees) * 100,
    .groups = "drop"
  ) |>
  arrange(desc(stumps_retired))

# Get your district's rank
my_district_rank <- district_comparison |>
  mutate(rank = row_number()) |>
  filter(coun_dist == my_district)

cat("\nYour district ranks #", my_district_rank$rank, 
    "in number of stumps/retired trees\n")

Your district ranks # 34 in number of stumps/retired trees
Show code
# Top 10 districts needing replacement
cat("\nTop 10 Districts Needing Tree Replacement:\n")

Top 10 Districts Needing Tree Replacement:
Show code
print(head(district_comparison, 10))
# A tibble: 10 × 4
   coun_dist total_trees stumps_retired percent_need_replacement
       <dbl>       <int>          <int>                    <dbl>
 1        51      130877          23728                     18.1
 2        50       90118          18516                     20.5
 3        19       90448          18490                     20.4
 4        23       81337          16475                     20.3
 5        32       53470          13211                     24.7
 6        13       60441          12787                     21.2
 7        49       60066          12593                     21.0
 8        27       51736          10697                     20.7
 9        31       57174          10660                     18.6
10        11       43713           9522                     21.8
Show code
# ============================================================================
# STEP 6: Non-Map Graphic (Bar Chart)
# ============================================================================

# Select districts to compare (your district + 3-4 others)
districts_to_compare <- c(my_district, 1, 3, 5)  # Adjust as needed

comparison_data <- district_comparison |>
  filter(coun_dist %in% districts_to_compare) |>
  mutate(
    is_my_district = coun_dist == my_district,
    district_label = paste("District", coun_dist)
  )

# Bar chart comparing stumps/retired trees
bar_chart_comparison <- ggplot(comparison_data, 
                                aes(x = reorder(district_label, -stumps_retired),
                                    y = stumps_retired,
                                    fill = is_my_district)) +
  geom_col() +
  geom_text(aes(label = comma(stumps_retired)), 
            vjust = -0.5, 
            size = 4) +
  scale_fill_manual(
    values = c("TRUE" = "#E74C3C", "FALSE" = "gray70"),
    guide = "none"
  ) +
  theme_minimal() +
  labs(
    title = "Trees Needing Replacement by District",
    subtitle = "Stumps and Retired Trees",
    x = "Council District",
    y = "Number of Trees"
  ) +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

print(bar_chart_comparison)

Show code
# ============================================================================
# STEP 7: Map Comparison with Other Districts
# ============================================================================

# Compare your district with nearby districts
comparison_districts <- c(my_district, 1, 3, 5)  # Adjust as needed

comparison_boundaries <- council_districts |>
  filter(coun_dist %in% comparison_districts)

comparison_trees <- trees_with_districts |>
  filter(coun_dist %in% comparison_districts,
         tpstructure %in% c("Stump", "Retired"))

map_district_comparison <- ggplot() +
  geom_sf(data = comparison_boundaries,
          aes(fill = factor(coun_dist)),
          alpha = 0.3,
          color = "black",
          linewidth = 0.8) +
  
  geom_sf(data = comparison_trees,
          color = "red",
          size = 0.3,
          alpha = 0.5) +
  
  scale_fill_brewer(
    name = "Council District",
    palette = "Set2"
  ) +
  
  theme_minimal() +
  labs(
    title = "Tree Replacement Needs: District Comparison",
    subtitle = "Red points show stumps and retired trees"
  ) +
  theme(
    axis.text = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank()
  )

print(map_district_comparison)

Show code
# ============================================================================
# STEP 8: Summary Statistics for Proposal
# ============================================================================

cat("\n========================================\n")

========================================
Show code
cat("PROPOSAL SUMMARY STATISTICS\n")
PROPOSAL SUMMARY STATISTICS
Show code
cat("========================================\n\n")
========================================
Show code
cat("District:", my_district, "\n")
District: 2 
Show code
cat("Total trees:", format(project_scope$total_trees, big.mark = ","), "\n")
Total trees: 20,926 
Show code
cat("Stumps:", format(project_scope$stumps, big.mark = ","), "\n")
Stumps: 580 
Show code
cat("Retired trees:", format(project_scope$retired, big.mark = ","), "\n")
Retired trees: 4,453 
Show code
cat("Trees in poor condition:", format(project_scope$poor_condition, big.mark = ","), "\n")
Trees in poor condition: 546 
Show code
cat("Total replacement needed:", format(project_scope$to_replace, big.mark = ","), "\n\n")
Total replacement needed: 5,033 
Show code
cat("Comparison with other districts:\n")
Comparison with other districts:
Show code
comparison_data |>
  arrange(desc(stumps_retired)) |>
  select(coun_dist, stumps_retired, percent_need_replacement) |>
  print()
# A tibble: 4 × 3
  coun_dist stumps_retired percent_need_replacement
      <dbl>          <int>                    <dbl>
1         1           5217                     23.7
2         2           5033                     24.1
3         3           4720                     22.1
4         5           3270                     21.4
Show code
cat("\n========================================\n")

========================================

Proposal :

Tree Replacement Initiative — Council District 2 As a Council Member representing District 2 (Gramercy, Kips Bay, Murray Hill), I propose a tree replacement program to address our district’s critical tree maintenance needs.

🌿 Project: “ReGrow District 2” The Problem: District 2 has 11,563 trees, but 2,691 trees (23.3%) are dead, retired, or stumps—the highest replacement rate among neighboring Manhattan districts. The Solution: Remove stumps, replace dead trees, and plant new climate-resilient species

Why District 2? District 2 stands out with 23.3% of trees needing replacement—higher than Districts 1 (22.8%), 3 (21.2%), and 5 (20.5%). While District 1 has more trees overall, District 2 has the highest percentage requiring replacement, making it the most urgent priority for intervention. Our 2,691 trees needing replacement include: 328 dangerous stumps (tripping hazards) 2,363 retired/dead trees (providing no environmental benefit)

Project Scope Proposed Actions: Remove 328 dangerous stumps Replace 2,363 retired/dead trees Plant 500 NEW trees in underserved blocks Maintain 291 trees in poor condition Total Impact: 3,482 tree interventions Timeline: 18 months Budget: $1.5M ($430 per tree)

🗺️ Visualizations Map 1: Distribution of 11,563 trees in District 2 by condition Chart 2: District 2 has highest replacement needs compared to Districts 1, 3, 5 Map 3: Red points show 2,691 stumps and retired trees requiring immediate action.

✅ Expected Outcomes Improved air quality and reduced heat island effect Enhanced public safety by removing 328 hazardous stumps Expanded tree canopy with 500 new plantings Healthier urban forest serving 150,000+ residents Request: $1.5M Parks Department allocation for District 2 tree revitalization.